TECHNICAL  REPORT  ARCCB-TR-95041 


AD 


LINEAR  TIME  ALGORITHM  FOR  POSITIVE  KERNEL 
SMOOTHING  WITH  APPLICATION  TO 
NONPARAMETRIC  PROBABILITY 
DENSITY  ESTIMATION 


ROYCE  W.  SOANES 


OCTOBER  1995 


US  ARMY  ARMAMENT  RESEARCH. 
DEVELOPMENT  AND  ENGINEERING  CENTER 

CLOSE  COMBAT  ARMAMENTS  CENTER 
BENET  LABORATORIES 
WATERVLIET,  N.Y.  12189-4060 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


1 9960320  034 


Dm;  qUALSTI  INSPECTED  I 


THIS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE 
COPY  FURNISHED  TO  DTIC 
CONTAINED  A  SIGNIFICANT 
NUMBER  OF  PAGES  WHICH  DO 
NOT  REPRODUCE  LEGIBLY. 


DISCLAIMER 


The  findings  in  this  report  ere  not  to  be  construed  as  an  official 
Department  of  the  Army  position  unless  so  designated  by  other  authorized 
documents . 

The  use  of  trade  name(s)  and/or  manufacturer(s)  does  not  constitute 
an  official  indorsement  or  approval. 


DESTRUCTION  NOTICE 

For  classified  documents,  follow  the  procedures  in  DoD  5200. 22-M, 
Industrial  Security  Manual,  Section  11*19  or  DoD  5200. 1-R,  Information 
Security  Program  Regulation,  Chapter  IX. 

For  unclassified,  limited  documents,  destroy  by  a.iy  method  that  will 
prevent  disclosure  of  contents  or  reconstruction  of  the  document. 

For  unclassified,  unlimited  documents,  destroy  when  the  report  is 
no  longer  needed.  Do  not  return  it  to  the  originator. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  I  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  anv  other  aspect  of  this 
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INTRODUCTION 


The  problem  we  are  concerned  with  in  this  report  is  developing  a  reliable  and  fast 
algorithm  for  continuous  smoothing  of  piecewise  linear  (or  piecewise  constant)  data  defined  over 
uniform  or  nonuniform  meshes.  The  theory  behind  what  we  will  do  is  covered  in  Reference  1. 
Here,  we  are  concerned  with  actual  computation.  Often,  large  amounts  of  empirical  data  are 
obtained  over  uniform  meshes  and  the  smoothing  is  carried  out  over  these  same  meshes.  This 
case  is  the  simplest  one  to  consider,  because  we  need  not  do  any  interpolation;  but  even  in  this 
simplest  of  cases,  we  must  guard  against  algorithmic  inefficiency  when  dealing  with  large  amounts 
of  data  consisting  of  perhaps  tens  of  thousands  of  points.  The  linear  time  property  of  an 
algorithm  is  therefore  quite  important  in  such  cases  and  enables  us  to  compute  in  seconds  what 
might  take  hours  otherwise.  The  situation  becomes  somewhat  more  complicated  when  the  data 
is  nonuniformly  spaced  and/or  when  we  wish  to  evaluate  the  smoothed  function  over  some  other 
mesh.  We  use  continuous  positive  kernel  smoothing  here,  because  of  its  shape  (monotonicity) 
preserving  properties  and  the  ease  with  which  we  may  interpolate  functional  and  derivative 
values. 


SMOOTHING  BY  REPEATED  AVERAGING 

In  this  section,  we  derive  the  smoothing  formulas  which  will  subsequently  be  applied. 
These  formulas  are  discussed  in  Reference  1,  but  we  include  their  derivation  here  for  the  sake  of 
convenience.  Consider  the  following  smoothing  operator  5: 

S  operating  on  function  /  at  point  x  is  simply  the  average  value  of  /  over  an  interval  of  length  2h 
with  x  as  the  center  of  the  interval.  The  parameter  h  is  called  the  window  parameter.  This 
integral  operator  can  also  be  written  in  kernel  form  as: 

if  \t\<h 
if  |f|;>A 


Applying  S  a  second  time  gives  us 


2 hJ*-k 


1  rx+h 
2 Wx-h 

1  Px+A  ft+h 


(2  hf 


i 


Changing  the  order  of  integration  in  this  double  integral  gives  us 

1  fx+2 h  rx*h. 


=  -^—[x  Ku)(.u+2h-x)du 
(2  h)2ixV> 

1  fx*Vi 


(2  hy 


fX+  fiju)(x+2h-u)du 


Writing  this  integral  in  kernel  form,  we  have 

S2Ax)  =  f  °J(u)K(u  -x)du 

J  —  CO 

if  -2h<u<0 

if0zu<2h 
{  0  if  \u\z.2h 

Note  that  the  discontinuous  kernel  corresponding  to  S  and  the  continuous  kernel  corresponding 
to  S2  are  both  nonnegative  or  essentially  positive  everywhere.  Although  we  will  not  prove  it  here, 
successive  applications  of  S  result  in  a  sequence  of  increasingly  smoother  and  wider  B-spline 
kernels  of  unit  area.  The  width  of  the  kernel  corresponding  to  Sk  is  2 kh.  We  now  leave  the 
subject  of  kernel  form  and  proceed  with  the  forms  we  actually  use  for  computation. 

Define  the  following  set  of  integrals  recursively: 

/;//»* 


*(«)= 


(2h+u) 


m2 

(2/t-tt) 

(2  hf 


where  a  is  an  arbitrary  value  in  the  domain  of  definition  of/.  Of  course,  the  derivative  of  any 
member  of  this  sequence  is  simply  the  preceding  member  of  the  sequence.  This  will  make  taking 
derivatives  a  trivial  operation. 
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As  before,  applying  S  to /once,  we  have 


=-^(/1(jc+A)-/1(jc-A)) 

2h 


Applying  5  a  second  time  and  noting  that  5  is  a  linear  operator,  we  have 


SVW- 


— (I2(x  +2h)  -2I2(x)  +I2(x  -2h)) 


Applying  5  a  third  time,  we  have 


Sftx) = — - —  (SI2(x +2h)  -2SIJx)  +SIJx-2h)) 
(2  hf 


-(/3(jc+3A)  -/,(* +A)  -2(/3(x+/z)  -/,(* -A))  +/,(x-A)  -L(x-3h)) 


-(IJx+3h)-3Ux+h)  +3L(x-h)  -Ux-3h)) 


We  can  then  derive  the  formula  for  the  fourth  smooth  or,  noting  the  appearance  of  the  binomial 
coefficients  with  alternating  signs,  we  can  simply  write  down  the  result  of  applying  S  a  fourth 
time: 

S4ftx)  =  — — r(J4(*  +4^)  _4/4(x  +2h)  +6/4(x)-4/4(x-2A)  +/4(x-4/t)) 

(2 h)* 
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Also  note  that  if  we  seek  the  derivative  of  this  fourfold  smoothed  function,  we  simply  reduce  the 
indices  of  the  Fs  by  unity  and  get 

—S*M=—(Ii(x+4h)-4Ux+2h)+6L(x)-4IJx-2h)+Ux-4h)) 
dx  (2  Kf 

These  ideas  are  elaborated  upon  in  Reference  1. 

COMPUTING  THE  INTEGRALS 

We  start  with  a  continuous  piecewise  linear  interpolation  of  the  data. 

Ay* 

/,(*) = y,+— (*-*,)  (Xt<.x<,  *j+1) 


For  Iv  we  have 


A(*)= f*M*= 

=/j (x)  +(X -x)y( + ^  -x)2  Ay‘ 


=  W  +yiAxi+-AxiAyi 


^(x^-Axfo+y^) 
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/2oo=  fxi1m=fXii1(t)dt+fxi1m 

Jx0  JxQ  Jxi 

+/x*  W 

=/2(x,)+/i(x4)Axi+Ax<?||yi+|y<+1j 


For  /j,  we  have 


^3C^+i)=^3C^)+/2^«)A*,+^/i<:^)A^f+^A^+^A>'«A^3 

»/^+^Ax,+|/l(x^Axf+A^-|yJ+^yl+1j 


This  is  as  far  in  the  sequence  of  Fs  as  we  will  go  in  this  report.  The  previous  set  of  formulas 
gives  us  two  things.  The  first  thing  is  the  recursion  relationships  we  must  use  to  compute  the 
nodal  values  of  the  first  three  integrals.  The  second  thing  is  the  formulas  for  interpolating  the 
integrals  at  arbitrary  points.  The  nodal  integral  values  are  essentially  computed  only  once  and 
the  amount  of  work  involved  is  obviously  proportional  to  the  amount  of  data  present.  The 
interpolations  may  obviously  be  done  in  constant  time,  once  the  intervals  of  the  various 
arguments  have  been  located.  Since,  in  the  evaluation  of  the  smoothed  function,  we  march 
through  the  data  only  once  as  we  locate  the  intervals  of  points  in  the  desired  mesh,  this 
component  of  the  algorithm  will  also  be  done  in  linear  time.  The  amount  of  work  done  to 
evaluate  the  smooth  approximation  will  therefore  be  kjn+k2m  arithmetic  operations,  where  n  is 
the  amount  of  data  and  m  is  the  number  of  points  at  which  we  wish  to  evaluate  the  smoothed 
function.  It  is  also  important  to  note  that  the  amount  of  work  involved  in  this  computational 
scheme  is  completely  independent  of  the  size  of  the  window  parameter  h. 


This  scheme  does  have  a  weakness,  however.  If  the  nodal  values  of  the  integrals  are 
computed  for  an  entire  very  large  set  of  data,  roundoff  error  in  the  integral  computation  and  in  the 
necessary  interpolations  (involving  subtractions)  will  completely  destroy  the  computation,  especially 
for  relatively  small  window  sizes.  The  results  will  be  so  contaminated  with  error  they  will 
meaningless.  For  a  small  to  intermediate  amount  of  data,  however,  the  results  will  be  very  good. 
It  is  therefore  wise  to  periodically  compute  the  nodal  values  of  the  integrals  in  smaller  clumps  as  we 
march  through  the  data.  Remember  that  the  lower  limit  of  integration  for  the  integrals  is 
analytically  arbitrary.  A  good  rule  of  thumb  is  to  compute  the  integrals  over  some  multiple  of  the 
total  window  width.  This  modification  naturally  forces  us  to  do  a  bit  of  extra  computation.  Each 
clump  of  integral  values  must  overlap  the  previous  clump  by  at  least  the  total  window  width. 

Therefore,  with  a  window  multiplier  of  ten,  we  will  do  only  about  ten  percent  more 
computational  work.  Sufficient  accuracy  can  still  be  attained  with  window  multipliers  as  large  as  one 
hundred.  In  any  case,  the  modified  algorithm  can  handle  any  amount  of  data  in  linear  time  with 
sufficient  accuracy.  An  involved  error  analysis  could  perhaps  be  done  for  this  algorithm,  but  it  is 
more  practical  to  evaluate  its  accuracy  by  applying  the  algorithm  to  an  artificial  set  of  data  obtained 
from  a  straight  line,  since  the  smoothing  operators  leave  such  data  undisturbed  analytically  but  not 
numerically.  All  we  need  do  is  to  compare  the  original  data  with  that  produced  by  the  algorithm 
with  its  inherent  roundoff  error.  Since  the  smoothing  formulas  require  data  on  either  side  of*  and 
the  domain  of  the  data  is  finite,  the  smoothing  operator  cannot  be  used  near  the  ends  of  the  data. 
For  this  case,  we  merely  use  the  simple  derivative  estimating  properties  of  the  smoothing  formulas 
to  expand  the  smoothed  function  in  a  Taylor  series  at  a  point  sufficiently  far  from  the  ends  and 
compute  the  smoothed  function  near  the  ends  via  these  Taylor  series.  Finally,  if  the  data  has  been 
taken  over  a  uniform  mesh  (equally  spaced  values  of  the  independent  variable)  and  the  smoothed 
function  is  desired  over  this  same  mesh,  only  the  nodal  values  of  the  integrals  are  necessary  and  no 
interpolations  need  be  done. 


NONPARAMETRIC  DENSITY  ESTIMATION 

The  simplest  and  most  common  example  of  nonparametric  density  estimation  is  the 
frequency  histogram.  Our  goal  in  this  section  is  to  improve  upon  the  discontinuous  histogram  by 
replacing  it  with  a  density  estimate  that  is  smooth  to  the  extent  of  being  twice  continuously 
differentiable.  Attainment  of  this  goal  will  give  us  a  density  that  is  not  only  more  aesthetically 
pleasing  to  the  eye,  but  one  with  which  we  can  make  more  accurate  probability  statements  and 
calculate  more  accurate  percentiles.  This  is  of  course  done  via  positive  kernel  smoothing.  Assuming 
that  a  sequence  of  values  of  a  random  variable  x  has  been  sorted  in  ascending  order,  we  can  define 
the  discrete  cumulative  distribution  function  as 


6 


w- 


0  ifx<xQ 


i+1  . 


ifxi<x<xM 

V 

1 


Note  that  this  function  is  defined  on  the  entire  real  line,  so  that  we  need  not  use  Taylor  series  to 
smooth  near  the  ends  of  a  finite  domain.  Also,  it  does  not  matter  how  the  distribution  function 
is  defined  at  the  data  points  or  whether  any  of  the  x  values  are  repeated  values.  Our  first 
estimate  of  the  density  will  be  the  derivative  of  the  fourth  smooth  of  Fn.  Hence,  we  need  three 
successive  integrals  of  Fn.  Since  the  values  of  x  are  always  nonuniformly  spaced,  interpolations 
are  mandatory. 

For  /j,  we  have 

[xFHm=pFnm+[x^dt 

Jxt>  *'xb  J*i  n 

=  I1(x)+—(x-xi) 

n 


(i+1)  Ax, 

For  I2,  we  have 

I2(x)=  f  x  I,m= f  Xh(x) + —{t-x)dt 

Jxq  Jxq  Jxi  n 

=  I2(x)  +Il(xi)(x-xt)  +-^(x-*i)2 
(i+l)Ax? 
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For  /3,  we  have 


i3(x)=  (xi2m=  (\m* fXcxj+ifyxt-xj+^it-xfdt 

Jxo  Jxi 

=  I}(x)  +I2(x)(x-xi)  +  ^-IfajXx  -x)2 + ^-(x  -x)3 
2  bn 

l  2  (i+l)Axf 

h(Xi+l)=  ^Xt^+^X^Xi+'^l^X^Xi  + 

Density  /  is  given  by 

Ax)  =  — —  U3(.x+4h)  -4I3(x+2h)  +6 13(x)  -4I3(x-2h)  +I3(x-4h)] 

(2 hf 

Note  that  although  we  have  derived  the  interpolation  formulas  for  /3  and  I2  in  passing,  only  the 
interpolation  formula  for  /3  is  necessary  since  we  have  no  need  of  Taylor  series  expansions.  Also, 
the  formula  for  I3(x)  can  be  used  to  extrapolate  beyond  the  last  point. 

PRESERVING  THE  VARIANCE 

The  sample  mean  computed  from  the  data  will  coincide  with  the  mean  computed  from 
the  density  estimate.  The  sample  variance  computed  from  the  data  will  not  coincide  with  the 
variance  computed  from  the  density  estimate.  The  density  variance  will  always  exceed  the 
sample  variance  by  an  ever  larger  amount  as  the  window  parameter  is  increased.  The  object  of 
this  section  is  to  derive  a  simple  transformation  which  will  keep  the  density  variance  the  same  as 
the  sample  variance  for  all  values  of  the  window  parameter  h.  We  now  return  to  the  use  of 
kernel  form  momentarily.  The  continuous  cumulative  distribution  function  estimate  is  given  by 

F(x)=SaFb(x)=  [yn(t)K(t-x)dt 

Integrating  by  parts,  we  have 

F(x)=  FnH udFn(i) 

1*1 

=  1--V  fXiK(u-x)du 

n (=0  _“ 

Differentiating  the  cumulative  distribution  function  gives  us  the  kernel  form  of  the  density: 

F'{x)=A> *)-i£  fXlK'(u-x)du=^£K(X-x) 

ft  i=o  “  ft  i=0 
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This  is  the  common  form  for  a  nonparametric  density  estimator  based  on  kernel  K.  However,  it  is 
exactly  the  wrong  form  to  use  for  computation.  If  it  were  used  to  evaluate  a  density  at  each  of  n 
points,  an  amount  of  work  proportional  to  the  square  of  the  amount  of  data  would  be  necessary. 
That  is,  instead  of  a  linear  time  algorithm,  we  would  have  a  quadratic  time  algorithm.  We  can, 
however,  use  this  formula  to  determine  the  relationship  between  the  variance  computed  from  the 
density,  the  sample  variance,  and  the  window  parameter  h.  Consider  the  expected  value  (via  f)  of 
an  arbitrary  function  p. 


E[p(x)]=  rPmx)dx=rp(x)*E  K{x-x)dx 

,»-i 

=~Y,  f“p(x+x.)K(x)dx 
ntZJ-K 


With  p(x)=x,  we  have 


£[*]=  V-f  =-£  f“(x+xf)K(x)dx=-^x.=ii 

ni=0J~°°  ni=t\ 


n-1 

nE 

ni= 0 


since  xK(x)  is  an  odd  function. 

We  therefore  see  that  the  sample  mean  and  the  density  computed  mean  are  the  same. 
Computing  the  variance  (again  via  /),  we  have 


m= oj=£[(x-p/  fjx+x-v/mdx 

--E  r(x2+2x(xr\if)+(xr\if)2)K(x)dx 
n£oJ-  1  1 

=/V*(x)d^£( xrnf 
J  —  nto  1 


1  i-0 

=S4x2(0)+d2 


But  SV  is  given  by 


S4*2 =— !-(/4Ct +4/z)  -4/4(x+2A) +6I4(x)  -4/4(x-2h)  +I4(x  -4h)) 
(2  Kf 


23-32-5 
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Hence 


S4x2(0) = — ^—~(2I4(4h) -8/4(2A)) 
(2  h)4 

1  [213*6  29,/t6  > 

(2A)4\23-32-5  23-32-5, 


=— (26-22)= 
32-5 


3 


We  therefore  finally  have  the  relation  between  the  density-computed  variance,  the  sample 
variance,  and  the  window  parameter: 


Now  let  g  be  the  density  of  a  new  random  variable  (*'). 

g(x)=  aj{\if+a(x-\if)) 


This  transformation  represents  a  simultaneous  compression  (about  the  mean)  off  in  the  x 
direction  and  a  stretching  of  /  in  the  y  direction.  Our  objective  now  is  to  compute  the 
transformation  parameter  a  in  terms  of  the  sample  variance  and  the  window  parameter  such  that 
the  variance  of  the  new  random  variable  computed  via  g  will  be  the  sample  variance. 

Considering  again  the  expected  value  (via  g)  of  an  arbitrary  function  p  of  the  new  random 
variable,  we  have 

E\p{x^\  =  f“p(x)g(x)dx 


=/~P(*)a/fy/+  a  (x-\i})dx 


=£V(*+ \i)ufi\if+ax)dx 


Letting  /?(*)= 1,  we  see  that  g  is  indeed  a  legitimate  density.  Letting />(*)=*,  we  see  that  the 
means  for/ and g  are  the  same.  Computing  the  variance  of*',  we  have 
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Therefore, 


p(x)=(x-\if)2 


Now,  since 

g(x)= 

or  equivalently 

g(”a^+kl/  =<X^X) 


we  compute  /  over  some  x  mesh  and  transform  the  density  data  according  to 

/-«/,  x-££+{x 
a 

in  order  to  obtain  the  g  density  data  that  preserves  both  the  sample  mean  and  the  sample  variance. 


OPTIMAL  WINDOW  PARAMETER 

In  this  section  we  develop  some  machinery  for  selecting  the  optimal  window  parameter  h  for 
a  nonparametric  density  estimator.  As  before,  we  call  the  density  estimator  g  and  we  will  call  the 
actual  density  y.  The  mean  square  error  (MSE)  for  a  given  value  of  x  is  given  by 

MSE(x)  =  £[(g(*)-y(;c))2]  =  £[(*(*) -£fe(x)]  +£fe(x)]  ~  Y(x))2] 

=E[(g(x)-E[gm2 

+2  (g(x)  -mxME[g(x)]-Y(x)) 

+(£fe(*)]-K*))2] 

=m(x)-mx)])2]<E[g(x)]-Y(x))2 

=V[g(x)]+B(x)2 

V\g(x)]  is  the  variance  ofg  and  B(x)  is  the  bias.  The  integral  of  the  mean  square  error  (IMSE)  is 
given  by 
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IMSE  =  f*~V[g(x)]+B(x)2dx 


Our  goal  is  to  approximately  minimize  IMSE  with  respect  to  the  window  parameter  h.  We  seek  a 
formula  that  is  asymptotically  valid  for  large  n  and  small  h.  Recalling  that 

g(x)= 


and 


n-l 

ni= 0 


where  K  is  the  kernel  resulting  from  four  applications  of  the  smoothing  operator  S,  we  have 

»-i 

gQc)=-YtK{xrlx-a{x-li)) 
n  <=o 


But 


and 


xi-fi-a(x~ii)=  xt-x+x-fi-a(x-fj,)=  x^x+Cl-aXx-fj) 


+—  (/r-0) 

3  o2 


We  therefore  have 


xrfi-a(x-fi )  -  xrx-h2P(x) 


where 


m= 


2(x-») 
3  o2 


Hence 

g(xy^£K(x-x-h2P) 

n  i=0 
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Since  the  x’s  are  identically  distributed  independent  random  variables,  the  expected  value  and 
variance  of  g  (via  y)  are  given  by 


mx)\=-'EEiK(xrx-h2py\ 

n  l= o 

=  aE[K(xrx-h2py\ 
ng(x)]=^£V[K(xrx-h2P)] 

n  i=0 

=  —VIKiXf-x-h2?)] 

n 


We  now  begin  to  obtain  simple  asymptotic  expressions  for  the  expectation  and  variance  of  these 
kernel  expressions.  Now  K  is  a  kernel  which  becomes  a  Dirac  delta  as  h  approaches  zero,  so  we 
take  note  of  the  fact  that 

hKQm>m, 

where  A:  is  a  standard  kernel  independent  of  h.  In  our  case,  the  support  of  k  is  (-4,4).  First,  we 
have 

E[K(xrx-h2p)]  =  [  *“K(u  -x  -h2P)  y(u)du 
= J  °°K(u)  y(u  +x+h  2 P)du  =  J  ~K(hu)  y(x +hu+h  2 P)hdu 
= f  *°°k(,u)Y(x+hu+h2  P)du 

Expanding  part  of  the  integrand  in  a  Taylor  series,  we  have 

Y(x +hu  +h  2  P)  =  y(x) + y\x){}w +h2P) 

+ —  Y//(x)(hu  +h  2  P)2 + —  y^XxXhu  +h  2  P)3 +0(h  4) 

2  6 

=  Y(x)+Y/(,x)(hu+h1p) 

+—Y//(x)(h2u2+2h3uP+h4p2) 

2 

+ — Y///(x)(h  3u 3  +3h  *u2p+3h  sup2 +h  6p3)  +0(h  4) 

6 

Since  k  is  an  even  function,  the  terms  involving  odd  powers  of  u  make  no  contribution  to  the 
integral.  The  contributing  terms  of  the  series  are  therefore 
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Hence, 


Therefore 


Yc(x+fui+h2P)=  y(x)  +h2fiY/(x)+—h2u2  y //(x)  +  0(h  4) 

2 


E[K(xrx-h2P)]=  f^m^Yix^p/ixyh^VXxyOih^ 

=  y<je>+  h2PY/(x)+—h2Y//(x)j*m>u2k(u)du+(Xh4) 

2  ** 

fywh 

■  i-jy^  -  \ 


E[K(xrx-h2P)]  =  y(x)+~(x-fi)  Y/(x)+^h2Y//(x)  +  0(A4) 
3  <r  3 

-  Y(x)+^-((x~n) y/(x)  +  o2x//(jc))+ £>(A4) 

3<r 


The  bias  is  therefore  given  approximately  by 


a  \y(x)+—((x-h)y/(x)+  o2y//(x))J-y(x) 

=(a-l  )y(x)+^-j-((x-h)y/(x)+  o2y//(x)) 

3  cr 

-  ^-z{y(x)Hx-i*)Ax)+  oVt*)) 

3a2 


In  order  to  obtain  the  variance  of  g,  we  must  first  compute 
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E^K^-x-h1  pf]  =  j  *K(u -x-h2P)2y(u)du  =  J  *K(u)2 y(u  +x +h 2fi)du 

=  f  *"°K(hu)2  y(x +hu  +h2fi)hdu  =  —  f“k(u)2Y(x+hu+h2  fi)du 
J  -tO  fcj  -M 

h  J—  h 


Now,  since 

V[K(xrx-h2p)\= E[K(x-x-h2P)2]-  E[K(x.-x-h2P)]2 
we  have  the  desired  variance  of  g 

'  *x)2+ 

cy(x)  yjx)2 
nh  n 

The  integral  of  the  mean  square  error  is  therefore 

IMSE  ~  4  f  *~Y(x)2dx+—  f  *“b(x)2dx 
nh  nJ-~  9<r¥J_“ 

Letting 

J  =J  “b(x)2dx 


and  setting  the  derivative  with  respect  to  h  of  the  integral  of  the  mean  square  error  equal  to 
zero,  we  have 


_ c  1 16ft3/ _q 

nh 2  9a 4 


Solving  for  h,  we  have 


Opt 


(  9co4\ 


i 

5 


[16nJj 


We  therefore  see  that  as  the  sample  size  becomes  large,  the  optimal  window  parameter  becomes 
small,  albeit  very  slowly.  We  now  simplify  the  integral  J.  First,  we  assume  that 

lim  xy(jc)=0=  lim  xy'ix) 

X-»±*»  X-±» 
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Hxf  =  {Y(.x)Hx-fi)Y/(x)+o2Y//(x)f 
= y(x)2 +(x-fi)2Y/(x)2+a4//(x)2 
+2 (x-n)y(x)Y/(x) 
+2o2y(x)y//(x) 
+2o2(x-h)y/(x)y//(x) 


Integrating  by  parts,  we  have 

Jt=f  \x - fi) y(x) Y/{x)dx  =^\x-y)Yix)dYix) 

=  -/  ~Y(x){Y(x)+(x-f*)Y/(x))dx  =  -f*~y(x)2dx-J \ 

and 

j2=fy*)Y"(x)dx  =fyxW(x)=  -fyixfdx 

and 


Ji=j*\x-n)Y/(x)Y//(x)dx  =  J+J(x -p)  Y/(x)d  y'(x) 

=  -/  *W(xt.Y/(x)  +(x-fi)  Y/\x))dx  =  -fy(x)2dx-J3 


from  which  we  get 


J=f  +“y(jt)2+(jt-/a)2  /(x)2 + a4y//(x)2d!*+2^-^  /J°K*)2<kj 

+2  a2 1- J  +Jy/W2di|+2a2|-|  / +  V(*)2<fr  j 

= J  +“((JC-m)2-3  a2)y/(x)2+  o4Y//(x)2dx 


We  can  get  a  clearer  idea  of  the  dependence  of  /  on  a  by  invoking  the  standardized  version  of  y. 
Let  T  be  the  standardized  version  with  zero  mean  and  unity  variance.  First,  we  have 


Hence, 


where  the  integral  factor  is  completely  independent  of  a.  We  therefore  have 


We  will  now  find  an  expression  for  optimal  h  for  a  particular  I\  A  candidate  for  this  particular 
case  must  naturally  be  twice  differentiable.  The  function  k  naturally  comes  to  mind  as  a  good 
unimodal  candidate,  having  zero  mean.  It  does  not,  however,  have  unity  variance  and  must 
therefore  be  properly  scaled.  Consider 

I\x)=ak(ax) 

f  °°x2IJx)dx=l  =  f  *“x2ak(qx)dx=^-  f  +“x2k(x)dx=-^- 
j-»  j-co  a2J~“  3  a2 


Hence, 


Also, 

r'(x)=  a2k'(ax),  r"(x)=  a3k"(ax) 


Computing  the  Jo  integral,  we  have 


Ja=j*“(x2-3)a4k/(ax)2+a6k//(ax)2dx 
3j  ~  (*!  -^k'(x)2+a2k"(x)2dx 
=—  ( *\9x2-36)k'(x)2+l6k//(x)2dx 


Therefore,  optimal  h  for  this  particular  standardized  density  is  given  by 

'  SU3[~Kx)2dx 


hcp,  =  a 


\1 

5 


32 n  f “(9x2-36)k'(x)2+16k"(xfdx 
\  Jo 


From  Reference  1,  we  have 


(4 -a;)3 -4(2 -s)3  jf  0s.,s2 


96 

if  2  *x  *4 
96 

0  if  x  >4 
k(-x)  if  x<0 


from  which  the  first  and  second  derivatives  of  k  can  be  obtained.  Then,  we  can  compute  the 
approximation 


This  asymptotic  approximation  can  be  used  as  a  nominal  estimate  of  optimal  h  for  a  unimodal 
density.  For  multimodal  data,  or  in  order  to  pick  out  more  of  the  variations  present  in  the  data, 
one  can  naturally  use  smaller  values  of  h,  but  we  at  least  have  a  value  to  start  with. 

Since  this  report  is  mainly  computation  oriented,  we  include  some  relevant  C  routines  in 
the  Appendix.  The  first  part  of  the  Appendix  is  devoted  to  the  C  preprocessor  commands  that 
describe  the  syntax  used  in  the  routines.  The  second,  third,  and  fourth  parts  of  the  Appendix  are 
devoted  to  uniform  mesh  smoothing,  nonparametric  density  estimation,  and  nonuniform  mesh 
smoothing,  respectively. 
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APPENDIX 
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#define  If  if( 

#define  Then  )  { 

#define  Elsif }  else  if( 

#define  Else  }  else  { 

#define  fl  }; 

#define  iF  ( 

#define  theN )  ? 

#define  elsE : 

#define  Fi ; 

#define  And  && 

#define  Or  1 1 
#def!ne  Do  for(;;)  { 

#define  oD  } 

#define  Un  if( 

#define  nU  )  break; 

#define  As  if(!( 

#define  sA ))  break; 

#define  NL  printf("\n"); 

#define  REAL(r)  printf(#r  "=%le  ”,r); 

#define  INT(i)  printf(#i  "=%ld  ”,i); 

#define  REALS(n,r)  {  long  _i_=0;  \ 
printf("\n");  \ 

Do  As  _i_  <  n  sA  \ 

printf(#r  ”[%ld]=%le  ”,_i_,r[_ij);  _L++;  oD  \ 
printf("\n");  }; 

#define  INTS(n,i)  {  long  J_=0;  \ 
printf("\n");  \ 

Do  As  J_  <  n  sA  \ 

printf(#i  ”[%ld]=%ld  ",  j  ,i[  ,i  ]);  J_++;  oD  \ 
printf("\n");  }; 

#define  AVER(bool,fname)  if(!(bool))  \ 

{  printf("\n"  #bool "  is  false  in  "  #fname);  \ 
exit(O);  }; 

#define  DENY(boolfname)  if(bool)  \ 

{  printf("\n"  #bool "  in  "  #fname);  \ 
exit(O);  }; 


#include  ”syntax.h" 

#define  DIM  500 

double  I1[DIM],I2[DIM],I3[DIM]; 


void  pkints(m,M,y) 
long  m,M; 
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double  *y; 

{ long  i,iml; 

i=m+l;  II [m] =0.0;  I2[m]=0.0;  I3[m]=0.0; 

Do  Un  i  >  M  nU 
iml=i-l; 

II  [i] =1 1  [im  1  ] + (y  [iml  ]  +y  [i])/2.0; 

I2[i] =I2[iml]+ (1 1  [im  1  ]+(y[iml]/3.0  +y[i]/6.0)); 
I3[i]=I3[iml]+(I2[iml]+(Il[iml]/2.0+(y[iml]/8.0+y[i]/24.0))); 
i++;  oD  } 


void  pksmooth(n,y,w,wm,s) 
long  n,w,wm; 
double  *y,*s; 

{  void  pkints(long, long, double*); 
double  H,f0,fl,f2,x; 

long  nint,m,M,i,imw,ipw,im3w,ip3w,j,ilast; 
DENY(n  >  DIM,pksmooth) 

DENY(6*w+l  >  n,pksmooth) 

DENY(wm  <  2,pksmooth) 
nint=(6*w+l)*wm; 

m=0;  M=m+nint;  If  M  >  =  n  Then  M=n-1;  fl 
pkints(m,M,y); 

H=2.0*(double)w;  H=H*H*H; 
i=3*w; 

iraw=i-w;  ipw=i+w; 

im3w=i-3*w;  ip3w=i+3*w; 

fD=(I3[ip3w]-I3[im3w]-3.0*(I3[ipw]-I3[imw]))/H; 

fl=(I2[ip3w]-I2[im3w]-3.0*(I2[ipw]-I2[imw]))/H; 

£2= (1 1  [ip3w]-1 1  [im3w]-3 .0*(1 1  [ipw]-Il  [imw]))/H; 

j=0; 

Do  Un  j  >  i  nU 
x=(double)(j-i); 
s[j]=f0+x*(fl+x*f2/2.0); 
j++;  oD 

i=j;  ilast=n-l-3*w; 

Do  imw=i-w;  ipw=i+w; 
im3w=i-3*w;  ip3w=i+3*w; 

If  ip3w  >  M 

Then  m=im3w;  M=m+nint; 

If  M  >  =  n  Then  M=n-1;  fl 
pkints(m,M,y);  fl 

fO = (13  [ip3w]-I3[im3w]-3.0*  (I3[ipw]-I3[imw]))/H; 
s[i]=fi); 

i++;  Un  i  >  ilast  nU  oD 
fl=(I2[ip3w]-I2[im3w]-3.0*(I2[ipw]-I2[imw]))/H; 
f2 = (1 1  [ip3w]-1 1  [im3w]-3 .0  *  (1 1  [ipw]-1 1  [imw]))/H; 
j=i;  i=ilast; 
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Do  As  j  <  n  sA 
x=(double)(j-i); 
s[j]=f0+x*(fl+x*£2/2.0); 
j++;  oD  } 


#include  "syntax.h" 

#include  "math.h” 

#define  DIM  500 

double  1 1  [DIM],I2[DIM],I3 [DIM] ; 

void  cdfints(n,x,m,M) 
long  n,m,M; 
double  *x; 

{  long  i,iml; 
double  dx,rn,ri; 

DENY(n  >  DIM,cdfints) 

If  m  <  0  Then  m=0;  fl 
If  M  >  =  n  Then  M=n-1;  fl 
DENY(m  >  =  M,cdfints) 

i=m+l;  ll[m]=0.0;  I2[m]=0.0;  I3[m]=0.0;  m=(double)n; 

Do  Un  i  >  M  nU 
iml=i-l;  dx=x[i]-x[iml]; 
ri=(double)i; 

Il[i]=Il[iml]+dx*ri/rn; 

I2[i]=I2[iml]+dx*(Il[iml]+dx*ri/(2.0*rn)); 

13  [i] = 13  [im  1  ] + dx*  (I2[im  1  ] + dx*  (1 1  [im  1  ]/2.0 + dx*  ri/(6.0*  m))); 
i++;  oD  } 


double  I3c(n,x,a,ila) 
long  n,ila; 
double  *x,a; 

{  double  dx,rn,rila; 

If  a  <  x[0]  Then  return  (0.0);  fl 
DENY(a  <  x[ila],I3c) 

dx=a-x[ila];  m=(double)n;  rila=(double)ila; 

return  (I3[ila]+dx*(I2[ila]+dx*(Il[ila]/2.0+dx*(rila+1.0)/(6.0*m)))); } 


void  npden(n,x,h,wm,nd,xd,yd) 
long  n,wm,nd; 
double  *x,h,*xd,*yd; 

{  void  cdfints(long, double*, long, long); 
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double  I3c(long, double  *  .double, long); 
void  musig(long,double*,double*, double*); 
double  xm,xl,xr,xll,xrr,xrrr,H,El,E2,mu,sig,alpha; 
long  id,im,il,ir,ill,irr,irrr,nml,i; 

DENY(n  >  DIM,npden) 

id=0;  im=0;  il=0;  ir=0;  ill=0;  irr=0;  nml=n-l; 

H=2.0*(double)h;  H=H*H*H*H; 
id=0;  irrr=0; 

Do  As  id  <  nd  sA 
xm=xd[id]; 

xl=xd[id]-2.0*h;  xr=xd[id]+2.0*h; 
xll  =xd[id]-4.0*h;  xrr =xd[id]  +4.0*h; 

Do  Un  im  ==  nml  nU  Un  x[im+l]  >  xm  nU  im++;  oD 
Do  Un  il  ==  nml  nU  Un  x[il+l]  >  xl  nU  il++;  oD 
Do  Un  ir  ==  nml  nU  Un  x[ir+l]  >  xr  nU  ir++;  oD 
Do  Un  ill  ==  nml  nU  Un  x[ill+l]  >  xll  nU  ill++;  oD 
Do  Un  irr  ==  nml  nU  Un  x[irr+l]  >  xrr  nU  irr++;  oD 
If  xrr  >  x[irrr]  And  irrr  <  nml 
Then  xrrr=xll+8.0*wm*h; 

Do  Un  irrr  ==  nml  nU  Un  xfirrr]  >  xrrr  nU  irrr++;  oD 
cdfints(n,x, ill, irrr);  fl 

El=I3c(n,x,xrr,irr)+6.0*I3c(n,x,xm,im)+I3c(n,x,xll,ill); 

E2=I3c(n,x,xr,ir)+I3c(n,x,xl,il); 

yd[id]=(El-4.0*E2)/H; 

If  yd[id]  <  0.0  Then  yd[id]=0.0;  fl 
id++;  oD 

musig(n,x,&mu,&sig); 

alpha=h/sig;  alpha=sqrt(1.0+4.0*alpha*alpha/3.0); 
i— 0; 

Do  As  i  <  nd  sA 
xd[i] = mu+ (xd[i]-mu)/alpha; 
yd[i]=alpha*yd[i]; 
i+  +;  oD  } 


void  musig(n,x,mu,sig) 
long  n; 

double  *x,*mu,*sig; 

{  long  i; 
double  s,d; 

i=0;  s=0.0;  Do  As  i  <  n  sA  s+=x[i];  i++;  oD 
*mu=s/(double)n; 

i=0;  s=0.0;  Do  As  i  <  n  sA  d=x[i]-(*mu);  s+=d*d;  i++;  oD 
*sig=sqrt(s/(double)n);  } 


void  sift(heap, root, last) 
double  *heap; 
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long  root, last; 

{  long  l,r; 
double  t; 

Do  l=2*root+l; 

If  1  >  last  Then  return;  fl 
r=l+l; 

If  r  >  last 

Then  If  heap[root]  <  heap[l] 

Then  t=heap[root];  heap[root]=heap[l];  heap[l]=t;  fl 
return;  fl 

If  heap[root]  >  =  heap[l]  And  heap[root]  >  =  heap[r] 

Then  return;  fl 
Ifheapfl]  >  heap[r] 

Then  t=heap[l];  heap[l]=heap[root];  heap[root]=t;  root=l; 

Else  t=heap[r];  heap[r]=heap[root];  heap[root]=t;  root=r;  fl  oD  } 


void  makeheap(array,last) 
double  *array; 
long  last; 

{  long  il,i2,i; 

void  sift(double*, long, long); 
il=0; 

Do  il=2*il+l;  Un  il  >  last  nU  oD 
il=(il-l)/2; 

Do  i2=il-l;  il=i2/2; 

If  last/2  <  i2  Then  i2=last/2;  fl 
i=il; 

Do  Un  i  >  i2  nU 

sift(array,i,last);  i++;  oD 
Un  il  ==  0  nU  oD  } 


void  heapsort(n,x) 
long  n; 
double  *x; 

{  long  last,i; 
double  t; 

void  makeheap(double*, long), sift(double*, long, long); 
last=n-l;  makeheap(x,last); 

Do  Un  last  ==  0  nU 
t=x[0];  x[0]=x[last];  x[last]=t;  last-; 
sift(x,0,last);  oD 
i=l; 

Do  As  i  <  n  sA 

DENY(x[i-l]  >  x[i],heapsort)  i++;  oD  } 
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#include  "syntax.h' 
#define  DIM  500 


double  II [DIM],I2[DIM],I3[DIM]; 

void  pkints(n,x,y,m,M) 
long  n,m,M; 
double  *x,*y; 

{ long  i,iml; 
double  dx; 

DENY(n  >  DIM,pkints) 

If  m  <  0  Then  m=0;  fl 

If  M  >  =  n  Then  M=n-1;  fl 

DENY(m  >=  M,pkints) 

i=m+l;  II [m] =0.0;  I2[m]=0.0;  I3[m]=0.0; 

Do  Un  i  >  M  nU 
iml=i-l;  dx=x[i]-x[iml]; 
Il[i]=Il[iml]+dx*(y[iml]+y[i])/2.0; 
I2[i]=I2[irnl]+dx*(Il[iml]+dx*(y[iml]/3.0+y[i]/6.0)); 
I3[i]=I3[iml]+dx*(I2[iml]+dx*(Il[iml]/2.0+dx* 
(y[iml]/8.0+y[i]/24.0))); 
i+  +;  oD  } 


double  Ilc(n,x,y,a,ila) 
long  n,ila; 
double  *x,*y,a; 

{  double  dx,r,E; 

If  ila  <  0  Then  ila=0;  fl 
If  ila  >  n-2  Then  ila=n-2;  fl 
dx=a-x[ila];  r=dx/(x[ila+l]-x[ila]); 
E= ((2.0-r)  *y[ila] + r  *y[ila + 1  ])/2.0; 
return  (Il[ila]+dx*E);  } 


double  I2c(n,x,y,a,ila) 
long  n,ila; 
double  *x,*y,a; 

{  double  dx,r,E; 

If  ila  <  0  Then  ila=0;  fl 
If  ila  >  n-2  Then  ila=n-2;  fl 
dx = a-x[ila];  r = dx/(x[ila + 1  ]-x[ila]); 
E=((3.0-r)*y[ila]+r*y[ila+l])/6.0; 
return  (I2[ila]+dx*(Il[ila]+dx*E)); } 
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double  I3c(n,x,y,a,ila) 
long  n,ila; 
double  *x,*y,a; 

{  double  dx,r,E; 

If  ila  <  0  Then  ila=0;  fl 

If  ila  >  n-2  Then  ila=n-2;  fl 

dx=a-x[ila];  r=dx/(x[ila+l]-x[ila]); 

E=((4.0-r)*y[ila]+r*y[ila+l])/24.0; 

return  (I3[ila]+dx*(I2[ila]+dx*(Il[ila]/2.0+dx*E))); } 


void  nupks(n,x,y,w,wm,ns,xs,ys) 
long  n,wm,ns; 
double  *x,*y,w,*xs,*ys; 

{  void  pkints(long, double*, double*, long, long); 
double  Ilc(long,double*, double*, double, long); 
double  I2c(long,double*, double*, double, long); 
double  I3c(long, double  *  ,doubl  e  *  .double, long); 
double  xl,xr,xll,xrr,xrrr,H,El ,E2,fD,fl,f2,xx; 
long  is,il,ir,ill,irr,irrr,j; 

DENY(n  >  DIM,nupks) 
is=0;  il=0;  ir=l;  ill=0;  irr=l; 

Do  Un  xs[is]-3.0*w  >=  x[0]  nU  is++;  oD 

xl=xs[is]-w;  xr=xs[is]+w;  xll=xs[is]-3.0*w;  xrr=xs[is]+3.0*w 

DENY(xrr  >  x[n-l],nupks) 

Do  Un  x[il+l]  >  xl  nU  il++;  oD 
Do  Un  x[ir]  >=  xr  nU  ir++;  oD 
Do  Un  x[ill+l]  >  xll  nU  ill++;  oD 
Do  Un  x[irr]  >=  xrr  nU  irr++;  oD 
irrr=irr;  xrrr=xll+6.0*wm*w; 

If  xrrr  >  x[n-l]  Then  xrrr=x[n-l];  fl 
Do  Un  x[irrr]  >=  xrrr  nU  irrr++;  oD 
pkints(n,x,y,ill,irrr); 

H =2.0* (double)w;  H=H*H*H; 
El=I3c(n,x,y,xrr,irr-l)-I3c(n,x,y,xll,ill); 
E2=I3c(n,x,y,xr,ir-l)-I3c(n,x,y,xl,il); 
fO=(El-3.0*E2)/H; 

El  =  I2c(n^c,y,xrr,irr-1  )-I2c(n,x,y,xll,ill) ; 

E2 = I2c(n,x,y,xr,ir- 1  )-I2c(n^,y,xl,il) ; 
fl =(E1-3.0*E2)/H; 

El=Ilc(n,x,yp{rr,irr-l)-Ilc(n^,y,xll,ili); 

E2=Ilc(n,x,y,xr,ir-l)-Ilc(n,x,y,xl,il); 

f2=(El-3.0*E2)/H; 

j=0; 

Do  As  j  <  is  sA 
xx=xs[j]-xs[is]; 
ys[j]=fD+xx*(fl +xx*f2/2.0); 
j++;  oD 
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Do  xl=xs[is]-w;  xr=xs[is]+w;  xll=xs[is]-3.0*w;  xrr=xs[is]+3.0*w; 
DENY (xrr  >  x[n-l],nupks) 

Do  Un  x[il+l]  >  xl  nU  il++;  oD 
Do  Un  x[ir]  >=  xr  nU  ir++;  oD 
Do  Un  x[ill+l]  >  xll  nU  ill++;  oD 
Do  Un  x[irr]  >  =  xrr  nU  irr++;  oD 
If  xrr  >  x[irrr] 

Then  xrrr=xll+6.0*wm*w; 

If  xrrr  >  x[n-l]  Then  xrrr=x[n-l];  fl 
Do  Un  xfirrr]  >=  xrrr  nU  irrr++;  oD 
pkints(n,x,y,ill,irrr);  fl 
El=I3c(n,x,y,xrr,irr-l)-I3c(n,x,y,xll,ill); 
E2=I3c(n,x,y,xr,ir-l)-I3c(n,x,y,xI,il); 
f0=(El-3.0*E2)/H; 
ys[is]=fO; 

Un  xs[is+l]+3.0*w  >  x[n-l]  nU  is++;  oD 
El=I2c(n,x,y,xrr,irr-l)-I2c(n,x,y,xll,in); 

E2=I2c(n,x,y,xr,ir-l)-I2c(n,x,y,xl,il); 

fl=(El-3.0*E2)/H; 

El=Ilc(n,x,y,xrr>irr-l)-Ilc(n,x,y,xll,ill); 

E2 = I  lc(n,x,y,xr,ir-l  )-1 1  c(n,x,y,xl,il); 

f2=(El-3.0*E2)/H; 

j=is+l; 

Do  As  j  <  ns  sA 
xx=xs[j]-xs[is]; 
ys[j]=fD+xx*(fl +xx*f2/2.0); 
j++;  oD } 
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